# Loading data
# experimental data
primary <- read.csv("../raw_data/primary_screen.csv")
speciesonly <- read.csv("../raw_data/species_only_data.csv")
indiv <- read.csv("../raw_data/individual_data.csv")
# re-create column for individualID
indiv <- indiv %>% # Create numbering variable
group_by(PaperID) %>%
mutate(IndividualID = paste(row_number(),unique(PaperID), sep="_"))
# host data
tree <- read.nexus("../raw_data/upham_tree_666.nex")
tax <- read.csv("../raw_data/taxonomy_mamPhy_5911species.csv")
bats <- tax[tax$ord=="CHIROPTERA",]
rodents <- tax[tax$ord=="RODENTIA",]
bat_tree <- drop.tip(tree, setdiff(tree$tip.label, bats$Species_Name))
rodent_tree <- drop.tip(tree, setdiff(tree$tip.label, rodents$Species_Name))
hosts <- read.csv("../raw_data/HostNames_NCBI_Upham.csv")
hosts$Host_Upham <- gsub(" ", "_", hosts$Host_Upham)
# virus names and taxonomy
vtax <- read.csv("../clean_data/Virus_taxonomy.csv")
viruses <- read.csv("../raw_data/VirusNames_translation_Feb23_2024.csv")
viruses <- left_join(viruses, vtax)
#subset to only host reported and upham
hosts <- hosts %>% select(Host_reported, Host_Upham) %>% unique()
# add names
primary <- dplyr::left_join(primary, viruses)
primary <- dplyr::left_join(primary, hosts)
speciesonly <- dplyr::left_join(speciesonly, viruses)
speciesonly <- dplyr::left_join(speciesonly, hosts)
indiv <- dplyr::left_join(indiv, viruses)
indiv <- dplyr::left_join(indiv, hosts)
# merge speciesonly with individual data
dat <- full_join(indiv, speciesonly)
# removing hosts and viruses with NA names
dat <- dat[!is.na(dat$Virus_ICTV),]
dat <- dat[!is.na(dat$Host_Upham),]
# make subset tree to sampled hosts
host_tree <- drop.tip(tree, setdiff(tree$tip.label, dat$Host_Upham))
# These are models of disease, so remove studies where pathology was not reported
dat <- dplyr::left_join(dat, unique(select(primary, PaperID, Virus_reported, Host_reported, PathologyReported_YN)))
# dim(dat)# 2923
dat <- dat[dat$PathologyReported_YN=="Y",]
# dim(dat)# 2726
# subset to only suscptible individuals / species
dat <- dat[dat$Susceptible_YN=="Y",]
# dim(dat) # 2238
# make virus tree
vtax <- as.data.frame(unclass(vtax), stringsAsFactors=TRUE)
frm <- ~superkingdom/realm/kingdom/phylum/class/order/family/genus/Virus_ICTV
vtree <- as.phylo(frm, data = vtax, collapse=FALSE)
vtree$edge.length <- rep(1, nrow(vtree$edge))
# include only viruses in subset data (e.g. after removing mole)
vtree <- drop.tip(vtree, setdiff(vtree$tip.label, dat$Virus_ICTV))
# plot(vtree)
# add host weight and cleaned dose and inoculation route data
dose_dat <- read.csv("../clean_data/dose_data.csv")
dat <- select(dat, -c(Dose_amount, Dose_unit))
dat <- dplyr::left_join(dat, dose_dat)
# adding host specificity data
host_range <- read.csv("../clean_data/host_range_data.csv")
dat <- dplyr::left_join(dat, host_range)
# dim(dat)# 2238
Visualization of brms default priors (blue) compared to custom priors (yellow) inspired by McElreath’s book Statistical Rethinking.
Priors for the slope parameters:
inv_logit <- function(x) {exp(x)/(1+exp(x))}
prior <- runif(n=10000, min=-100, max=100) # brms default is actually -Inf to +Inf
p <- inv_logit(prior)
prior2 <- rnorm(n=10000, 0, 1.5) # McElreath
p2 <- inv_logit(prior2)
plot(density(p, adj=0.1), xlim=c(0,1), col="#56B4E9", main="", xlab=expression(paste("inv_logit (",beta,")")))
lines(density(p2, adj=0.1), col="#E69F00")
The default improper uniform prior used by brms is flat on the log-odds scale, so when converted to the probabilty scale, puts excess mass near 0 and 1. A Normal(0, 1.5) prior is more reasonable, having a flatter distribution with a mode at 0.5 on the probablity scale.
Priors for the family-level variance parameters:
prior <- rstudent_t(n=10000, df=3, mu=0, sigma=10) # brms default
prior <- prior[prior>0]
# p <- inv_logit(prior)
prior2 <- rnorm(n=10000, 0, 1) # McElreath
prior2 <- prior2[prior2>0]
# p2 <- inv_logit(prior2)
plot(density(prior2, adj=0.1), xlim=c(0,10), col="#E69F00", main="", xlab=expression(paste(sigma)))
lines(density(prior, adj=0.1), col="#56B4E9")
The brms prior is extremely long tailed. Considering we are fitting a non-linear model in which changes in log-odds over over 4 have a diminishing influence due to the ceiling effect of the binomial distribution. With this in mind, a prior of Normal(0,1) constrains the posterior to a more realistic range, and is likely to improve sampling efficiency.
# create a host-virus-study level dataset
spec <- dat
spec <- spec[spec$Susceptible_YN==1,]
spec <- spec[!is.na(spec$Susceptible_YN),]
spec %>% group_by(PaperID, Virus_ICTV, Host_Upham, Host_order, Reservoir_match, ses_pd, EvoIso) %>% summarise(n_indivID=n_distinct(IndividualID), N_individuals=sum(unique(N_individuals), na.rm=TRUE)) %>% View()
# any disease, maximum severity, any mortality, number of susceptible individuals
spec <- spec %>% group_by(PaperID, Virus_ICTV, Host_Upham, Host_order, Reservoir_match, ses_pd, EvoIso) %>%
summarise(n_individuals = sum(
c(n_distinct(IndividualID, na.rm=TRUE), unique(N_individuals)), na.rm=TRUE),
n_studies = n_distinct(PaperID),
Disease_YN=max(Disease_YN, na.rm=TRUE),
severity_rank=max(severity_rank, na.rm=TRUE),
Mortality_YN=max(Mortality_YN, na.rm=TRUE)) %>% unique()
# for combinations in which mortality was all NA max() returns -Inf
# change this back to NA
spec$Mortality_YN[is.infinite(spec$Mortality_YN)] <- NA
# table(spec$Disease_YN)
# 123/(70+123) # 63.7% of h-v-s combinations reported disease
range(spec$n_individuals, na.rm=TRUE)
# 1 - 492
# additional columns for modelling separate phylo and non-phylo species-level effects
spec$Host_name <- spec$Host_Upham
spec$Virus_name <- spec$Virus_ICTV
# scaling / transform continuous predictors
spec$ses_pd <- scale_half(spec$ses_pd)
spec$EvoIso <- scale_half(spec$EvoIso)
spec$n_individuals <- scale_half(log(spec$n_individuals))
if (!file.exists("../fit_models/spec_disease_base.rds")) {
spec_disease_m0 <- brm(Disease_YN ~ Host_order,
data=spec, family=bernoulli(link = "logit"),
prior=c(prior(normal(0, 1), class = Intercept),
prior(normal(0, 1.5), class = b)),
iter=2000, thin=1,
control=list(adapt_delta=0.8, max_treedepth=10), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_disease_m0, "../fit_models/spec_disease_base.rds")
} else { spec_disease_m0 <- readRDS("../fit_models/spec_disease_base.rds")}
model <- spec_disease_m0
# summary(model, prob=0.9)
# save memory
# rm(spec_disease_m0)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Mean | CI | CI_low | CI_high | Rhat |
|---|---|---|---|---|---|
| b_Intercept | 0.07 | 0.90 | -0.25 | 0.40 | 1.00 |
| b_Host_orderChiroptera | 1.02 | 0.90 | 0.53 | 1.52 | 1.00 |
p1 <- cond_eff_plot(model)
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
Sample size: 193
With no hierarchical effects and only the bat/rodent predictor, the model suggests that bats are more likely than rodents to suffer overt disease.
However, once include hierarchical effects for host and virus species, this effect goes away:
if (!file.exists("../fit_models/spec_disease_base_phylo.rds")) {
spec_disease_m1 <- brm(Disease_YN ~ Host_order +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=3000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.8, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_disease_m1, "../fit_models/spec_disease_base_phylo.rds")
} else { spec_disease_m1 <- readRDS("../fit_models/spec_disease_base_phylo.rds")}
model <- spec_disease_m1
# summary(model, prob=0.9)
# save memory
# rm(spec_disease_m1)
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.63 | 0.90 | -0.68 | 1.94 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.06 | 0.90 | -1.44 | 1.30 | 1.00 | |
| sd_Host_name__Intercept | random | 0.56 | 0.90 | 0.05 | 1.32 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.81 | 0.90 | 0.08 | 1.78 | 1.00 | Host_Upham |
| sd_Virus_ICTV__Intercept | random | 1.44 | 0.90 | 0.30 | 2.58 | 1.01 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.66 | 0.90 | 0.64 | 2.60 | 1.01 | Virus_name |
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-10,10)) +
ylab("Conditional log-odds of disease") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
ggsave("../plots_tables/base_phylo_species_level_disease_comboplot.pdf", width=8, height=3)
Sample size: 193
if (!file.exists("../fit_models/spec_disease_full.rds")) {
spec_disease_full <-
brm(Disease_YN ~ Host_order + n_individuals +
ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.92, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_disease_full, "../fit_models/spec_disease_full.rds")
} else { spec_disease_full <- readRDS("../fit_models/spec_disease_full.rds")}
model <- spec_disease_full
# summary(model, prob=0.9)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.98 | 0.90 | -0.35 | 2.29 | 1.00 | |
| b_Host_orderChiroptera | fixed | -9.55e-03 | 0.90 | -1.54 | 1.50 | 1.00 | |
| b_n_individuals | fixed | 1.27 | 0.90 | 0.06 | 2.47 | 1.00 | |
| b_ses_pd | fixed | -0.69 | 0.90 | -2.17 | 0.82 | 1.00 | |
| b_EvoIso | fixed | 0.39 | 0.90 | -0.71 | 1.50 | 1.00 | |
| sd_Host_name__Intercept | random | 0.57 | 0.90 | 0.05 | 1.37 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.80 | 0.90 | 0.07 | 1.82 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 1.74 | 0.90 | 0.92 | 2.64 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 1.17 | 0.90 | 0.12 | 2.40 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.14 | 0.90 | 0.15 | 2.23 | 1.00 | Virus_name |
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) +
ylab("Conditional log-odds of disease") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_disease_full_comboplot.pdf", width=8, height=3)
Sample size: 193
if (!file.exists("../fit_models/spec_mort_full.rds")) {
spec_mort_full <-
brm(Mortality_YN ~ Host_order + n_individuals +
ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_mort_full, "../fit_models/spec_mort_full.rds")
} else { spec_mort_full <- readRDS("../fit_models/spec_mort_full.rds")}
model <- spec_mort_full
# summary(model, prob=0.9)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.05 | 0.90 | -1.13 | 1.20 | 1.00 | |
| b_Host_orderChiroptera | fixed | 0.49 | 0.90 | -1.05 | 2.07 | 1.00 | |
| b_n_individuals | fixed | 2.34 | 0.90 | 1.00 | 3.71 | 1.00 | |
| b_ses_pd | fixed | -1.64 | 0.90 | -2.92 | -0.35 | 1.00 | |
| b_EvoIso | fixed | 1.03 | 0.90 | -0.05 | 2.15 | 1.00 | |
| sd_Host_name__Intercept | random | 0.75 | 0.90 | 0.07 | 1.66 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.93 | 0.90 | 0.12 | 1.93 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 1.27 | 0.90 | 0.33 | 2.19 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 0.71 | 0.90 | 0.06 | 1.79 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.79 | 0.90 | 0.07 | 1.78 | 1.00 | Virus_name |
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) +
ylab("Conditional log-odds of mortality") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar2 <- intervals_plot(model)
plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_mort_full_comboplot.pdf", width=8, height=3)
Sample size: 154
if (!file.exists("../fit_models/spec_severity_full.rds")) {
spec_severity_full <-
brm(severity_rank ~ Host_order + n_individuals +
ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec, family=cumulative(probit),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.95, max_treedepth=14), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_severity_full, "../fit_models/spec_severity_full.rds")
} else { spec_severity_full <- readRDS("../fit_models/spec_severity_full.rds")}
model <- spec_severity_full
# summary(model, prob=0.9)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept[1] | fixed | -0.97 | 0.90 | -1.88 | -0.09 | 1.00 | |
| b_Intercept[2] | fixed | -0.53 | 0.90 | -1.42 | 0.35 | 1.00 | |
| b_Intercept[3] | fixed | -0.25 | 0.90 | -1.13 | 0.64 | 1.00 | |
| b_Intercept[4] | fixed | 0.49 | 0.90 | -0.40 | 1.41 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.25 | 0.90 | -1.49 | 0.97 | 1.00 | |
| b_n_individuals | fixed | 1.46 | 0.90 | 0.74 | 2.27 | 1.00 | |
| b_ses_pd | fixed | -0.66 | 0.90 | -1.65 | 0.32 | 1.00 | |
| b_EvoIso | fixed | 0.32 | 0.90 | -0.35 | 1.00 | 1.00 | |
| sd_Host_name__Intercept | random | 0.57 | 0.90 | 0.06 | 1.29 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.80 | 0.90 | 0.09 | 1.64 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 0.94 | 0.90 | 0.51 | 1.43 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 0.95 | 0.90 | 0.16 | 1.86 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.77 | 0.90 | 0.10 | 1.47 | 1.00 | Virus_name |
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) +
ylab("Conditional log-odds of severity") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar3 <- intervals_plot(model)
plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_severity_full_comboplot.pdf", width=8, height=3)
Sample size: 193
## Joint plot
global_ylim <- scale_y_continuous(lim=c(-12,12))
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))
joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)
ggsave("../plots_tables/species_level_conditionals_intervals_full.pdf", joint_plot, width=10.5, height=5.5)
Sample size: 193
if (!file.exists("../fit_models/spec_disease_resMatch.rds")) {
spec_disease_resMatch <- brm(Disease_YN ~ Host_order + n_individuals +
Reservoir_match + ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_disease_resMatch, "../fit_models/spec_disease_resMatch.rds")
} else { spec_disease_resMatch <- readRDS("../fit_models/spec_disease_resMatch.rds")}
model <- spec_disease_resMatch
# summary(model, prob=0.9)
# save memory
rm(spec_disease_resMatch)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.98 | 0.90 | -0.59 | 2.56 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.41 | 0.90 | -1.96 | 1.12 | 1.00 | |
| b_n_individuals | fixed | 1.02 | 0.90 | -0.12 | 2.21 | 1.00 | |
| b_Reservoir_match | fixed | 0.09 | 0.90 | -1.30 | 1.45 | 1.00 | |
| b_ses_pd | fixed | -0.27 | 0.90 | -1.83 | 1.30 | 1.00 | |
| b_EvoIso | fixed | 0.26 | 0.90 | -0.92 | 1.42 | 1.00 | |
| sd_Host_name__Intercept | random | 0.57 | 0.90 | 0.04 | 1.40 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.80 | 0.90 | 0.07 | 1.82 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 1.44 | 0.90 | 0.47 | 2.42 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 1.50 | 0.90 | 0.21 | 2.82 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.22 | 0.90 | 0.16 | 2.36 | 1.00 | Virus_name |
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) +
ylab("Conditional log-odds of disease") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_disease_resMatch_comboplot.pdf", width=8, height=3)
Sample size: 176
if (!file.exists("../fit_models/spec_mort_resMatch.rds")) {
spec_mort_resMatch <- brm(Mortality_YN ~ Host_order + n_individuals +
Reservoir_match + ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_mort_resMatch, "../fit_models/spec_mort_resMatch.rds")
} else { spec_mort_resMatch <- readRDS("../fit_models/spec_mort_resMatch.rds")}
model <- spec_mort_resMatch
# summary(model, prob=0.90)
# save memory
rm(spec_mort_resMatch)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Mortality_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.12 | 0.90 | -1.28 | 1.48 | 1.00 | |
| b_Host_orderChiroptera | fixed | 0.47 | 0.90 | -1.10 | 2.03 | 1.00 | |
| b_n_individuals | fixed | 2.10 | 0.90 | 0.82 | 3.42 | 1.00 | |
| b_Reservoir_match | fixed | -0.20 | 0.90 | -1.52 | 1.12 | 1.00 | |
| b_ses_pd | fixed | -1.38 | 0.90 | -2.82 | 0.06 | 1.00 | |
| b_EvoIso | fixed | 0.94 | 0.90 | -0.25 | 2.13 | 1.00 | |
| sd_Host_name__Intercept | random | 0.75 | 0.90 | 0.06 | 1.71 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.93 | 0.90 | 0.10 | 1.94 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 1.20 | 0.90 | 0.26 | 2.18 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 0.80 | 0.90 | 0.07 | 1.95 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.93 | 0.90 | 0.10 | 2.02 | 1.00 | Virus_name |
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) +
ylab("Conditional log-odds of mortality") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar2 <- intervals_plot(model)
plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_mortality_resMatch_comboplot.pdf", width=8, height=3)
Sample size: 139
if (!file.exists("../fit_models/spec_severity_resMatch.rds")) {
spec_severity_resMatch <- brm(severity_rank ~ Host_order + n_individuals +
Reservoir_match + ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec, family=cumulative(probit),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.95, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_severity_resMatch, "../fit_models/spec_severity_resMatch.rds")
} else { spec_severity_resMatch <- readRDS("../fit_models/spec_severity_resMatch.rds")}
model <- spec_severity_resMatch
# summary(model, prob=0.90)
# save memory
rm(spec_severity_resMatch)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$severity_rank, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept[1] | fixed | -1.02 | 0.90 | -2.10 | 0.03 | 1.00 | |
| b_Intercept[2] | fixed | -0.67 | 0.90 | -1.75 | 0.37 | 1.00 | |
| b_Intercept[3] | fixed | -0.36 | 0.90 | -1.44 | 0.70 | 1.00 | |
| b_Intercept[4] | fixed | 0.39 | 0.90 | -0.68 | 1.48 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.39 | 0.90 | -1.63 | 0.87 | 1.00 | |
| b_n_individuals | fixed | 1.29 | 0.90 | 0.58 | 2.09 | 1.00 | |
| b_Reservoir_match | fixed | -0.10 | 0.90 | -1.04 | 0.81 | 1.00 | |
| b_ses_pd | fixed | -0.35 | 0.90 | -1.47 | 0.79 | 1.00 | |
| b_EvoIso | fixed | 0.20 | 0.90 | -0.60 | 0.99 | 1.00 | |
| sd_Host_name__Intercept | random | 0.63 | 0.90 | 0.05 | 1.39 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.76 | 0.90 | 0.09 | 1.59 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 0.89 | 0.90 | 0.43 | 1.39 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 1.15 | 0.90 | 0.19 | 2.15 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.90 | 0.90 | 0.15 | 1.72 | 1.00 | Virus_name |
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) +
ylab("Conditional log-odds of severity") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar3 <- intervals_plot(model)
plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_severity_resMatch_comboplot.pdf", width=8, height=3)
Sample size: 176
## Joint plot
global_ylim <- scale_y_continuous(lim=c(-12,12))
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))
joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)
ggsave("../plots_tables/species_level_conditionals_intervals_resMatch.pdf", joint_plot, width=10.5, height=5.5)
Sample size: 176
lyssaviruses <- viruses$Virus_ICTV[viruses$genus=="Lyssavirus"] %>% unique()
spec_nolys <- spec[!spec$Virus_ICTV%in%lyssaviruses,]
if (!file.exists("../fit_models/spec_disease_nolys.rds")) {
spec_disease_nolys <- brm(Disease_YN ~ Host_order + n_individuals +
ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec_nolys, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_disease_nolys, "../fit_models/spec_disease_nolys.rds")
} else { spec_disease_nolys <- readRDS("../fit_models/spec_disease_nolys.rds")}
model <- spec_disease_nolys
# summary(model, prob=0.9)
# # save memory
rm(spec_disease_nolys)
# # pp <- brms::posterior_predict(model)
# # ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
# param_tab(model)
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditonal log-odds of disease") +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_disease_nolys.pdf", width=8, height=3)
Sample size: 148
if (!file.exists("../fit_models/spec_mort_nolys.rds")) {
spec_mort_nolys <- brm(Mortality_YN ~ Host_order + n_individuals +
ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec_nolys, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_mort_nolys, "../fit_models/spec_mort_nolys.rds")
} else { spec_mort_nolys <- readRDS("../fit_models/spec_mort_nolys.rds")}
model <- spec_mort_nolys
# summary(model, prob=0.9)
# save memory
rm(spec_mort_nolys)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Mortality_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.14 | 0.90 | -0.97 | 1.26 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.10 | 0.90 | -1.69 | 1.53 | 1.00 | |
| b_n_individuals | fixed | 1.94 | 0.90 | 0.49 | 3.43 | 1.00 | |
| b_ses_pd | fixed | -1.62 | 0.90 | -3.18 | -0.09 | 1.00 | |
| b_EvoIso | fixed | 1.33 | 0.90 | 0.16 | 2.55 | 1.00 | |
| sd_Host_name__Intercept | random | 0.66 | 0.90 | 0.05 | 1.55 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.75 | 0.90 | 0.07 | 1.73 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 1.39 | 0.90 | 0.46 | 2.34 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 0.67 | 0.90 | 0.05 | 1.71 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.76 | 0.90 | 0.06 | 1.76 | 1.00 | Virus_name |
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditional log-odds of mortality") +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar2 <- intervals_plot(model)
plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_mortality_nolys_comboplot.pdf", width=8, height=3)
Sample size: 113
if (!file.exists("../fit_models/spec_severity_nolys.rds")) {
spec_severity_nolys <- brm(severity_rank ~ Host_order + n_individuals +
ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec_nolys, family=cumulative(probit),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_severity_nolys, "../fit_models/spec_severity_nolys.rds")
} else { spec_severity_nolys <- readRDS("../fit_models/spec_severity_nolys.rds")}
model <- spec_severity_nolys
# summary(model, prob=0.9)
# save memory
rm(spec_severity_nolys)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$severity_rank, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept[1] | fixed | -1.06 | 0.90 | -1.91 | -0.23 | 1.00 | |
| b_Intercept[2] | fixed | -0.58 | 0.90 | -1.40 | 0.24 | 1.00 | |
| b_Intercept[3] | fixed | -0.29 | 0.90 | -1.10 | 0.52 | 1.00 | |
| b_Intercept[4] | fixed | 0.35 | 0.90 | -0.46 | 1.19 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.41 | 0.90 | -1.66 | 0.93 | 1.00 | |
| b_n_individuals | fixed | 1.23 | 0.90 | 0.41 | 2.07 | 1.00 | |
| b_ses_pd | fixed | -1.04 | 0.90 | -2.15 | 0.02 | 1.00 | |
| b_EvoIso | fixed | 0.46 | 0.90 | -0.20 | 1.14 | 1.00 | |
| sd_Host_name__Intercept | random | 0.56 | 0.90 | 0.05 | 1.34 | 1.01 | Host_name |
| sd_Host_Upham__Intercept | random | 0.85 | 0.90 | 0.10 | 1.75 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 0.94 | 0.90 | 0.49 | 1.44 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 0.45 | 0.90 | 0.04 | 1.16 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.84 | 0.90 | 0.16 | 1.53 | 1.00 | Virus_name |
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditional log-odds of severity") +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar3 <- intervals_plot(model)
plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_severity_nolys_comboplot.pdf", width=8, height=3)
Sample size: 148
## Joint plot
global_ylim <- scale_y_continuous(lim=c(-12,12))
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))
joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)
ggsave("../plots_tables/species_level_conditionals_intervals_nolys.pdf", joint_plot, width=10.5, height=5.5)
spec_heteros <- spec[spec$Reservoir_match%in%0,]
if (!file.exists("../fit_models/spec_disease_heteroHosts.rds")) {
spec_disease_heteros <- brm(Disease_YN ~ Host_order + n_individuals +
ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec_heteros, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=10), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_disease_heteros, "../fit_models/spec_disease_heteroHosts.rds")
} else { spec_disease_heteros <- readRDS("../fit_models/spec_disease_heteroHosts.rds")}
model <- spec_disease_heteros
# summary(model, prob=0.9)
# # save memory
rm(spec_disease_heteros)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
# param_tab(model)
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditonal log-odds of disease") +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(ar1)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_disease_heteroHosts_comboplot.pdf", width=8, height=3)
Sample size: 81
if (!file.exists("../fit_models/spec_mort_heteroHosts.rds")) {
spec_mort_heteros <- brm(Mortality_YN ~ Host_order + n_individuals +
ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec_heteros, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=10), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_mort_heteros, "../fit_models/spec_mort_heteroHosts.rds")
} else { spec_mort_heteros <- readRDS("../fit_models/spec_mort_heteroHosts.rds")}
model <- spec_mort_heteros
# summary(model, prob=0.9)
# save memory
rm(spec_mort_heteros)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Mortality_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.40 | 0.90 | -0.87 | 1.66 | 1.00 | |
| b_Host_orderChiroptera | fixed | 0.16 | 0.90 | -1.51 | 1.84 | 1.00 | |
| b_n_individuals | fixed | 2.05 | 0.90 | 0.35 | 3.74 | 1.00 | |
| b_ses_pd | fixed | -0.91 | 0.90 | -2.78 | 0.94 | 1.00 | |
| b_EvoIso | fixed | 0.23 | 0.90 | -1.28 | 1.74 | 1.00 | |
| sd_Host_name__Intercept | random | 0.92 | 0.90 | 0.08 | 2.12 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.69 | 0.90 | 0.06 | 1.70 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 0.90 | 0.90 | 0.08 | 2.06 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 0.64 | 0.90 | 0.05 | 1.64 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.91 | 0.90 | 0.08 | 2.03 | 1.00 | Virus_name |
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditional log-odds of mortality") +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar2 <- intervals_plot(ar2)
plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_mortality_heteroHosts_comboplot.pdf", width=8, height=3)
Sample size: 61
if (!file.exists("../fit_models/spec_severity_heteroHosts.rds")) {
spec_severity_heteros <- brm(severity_rank ~ Host_order + n_individuals +
ses_pd + EvoIso +
(1|PaperID) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=spec_heteros, family=cumulative(probit),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=10), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(spec_severity_heteros, "../fit_models/spec_severity_heteroHosts.rds")
} else { spec_severity_heteros <- readRDS("../fit_models/spec_severity_heteroHosts.rds")}
# pairs(spec_severity_heteros)
model <- spec_severity_heteros
# summary(model, prob=0.9)
# save memory
rm(spec_severity_heteros)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$severity_rank, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept[1] | fixed | -1.15 | 0.90 | -2.18 | -0.15 | 1.00 | |
| b_Intercept[2] | fixed | -1.00 | 0.90 | -2.03 | -7.83e-03 | 1.00 | |
| b_Intercept[3] | fixed | -0.57 | 0.90 | -1.57 | 0.40 | 1.00 | |
| b_Intercept[4] | fixed | 0.23 | 0.90 | -0.77 | 1.24 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.40 | 0.90 | -1.74 | 0.97 | 1.00 | |
| b_n_individuals | fixed | 1.67 | 0.90 | 0.44 | 3.01 | 1.00 | |
| b_ses_pd | fixed | -1.10 | 0.90 | -2.77 | 0.54 | 1.00 | |
| b_EvoIso | fixed | 0.05 | 0.90 | -1.20 | 1.35 | 1.00 | |
| sd_Host_name__Intercept | random | 1.42 | 0.90 | 0.47 | 2.39 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.55 | 0.90 | 0.04 | 1.40 | 1.00 | Host_Upham |
| sd_PaperID__Intercept | random | 0.58 | 0.90 | 0.05 | 1.35 | 1.00 | PaperID |
| sd_Virus_ICTV__Intercept | random | 0.70 | 0.90 | 0.05 | 1.75 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.34 | 0.90 | 0.43 | 2.30 | 1.00 | Virus_name |
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditional log-odds of severity") +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar3 <- intervals_plot(model)
plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/species_level_severity_heteroHosts_comboplot.pdf", width=8, height=3)
Sample size: 81
## Joint plot
global_ylim <- scale_y_continuous(lim=c(-12,12))
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))
joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)
ggsave("../plots_tables/species_level_conditionals_intervals_heteroHosts.pdf", joint_plot, width=10.5, height=5.5)
Sample size: 81
| Parameter | Effects | Component | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | conditional | 1.22 | 0.90 | -0.64 | 3.09 | 1.00 | |
| b_Host_orderChiroptera | fixed | conditional | 0.64 | 0.90 | -1.30 | 2.43 | 1.00 | |
| sd_Dose_unit__Intercept | random | conditional | 1.98 | 0.90 | 1.11 | 3.02 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | conditional | 2.05 | 0.90 | 1.52 | 2.58 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | conditional | 1.87 | 0.90 | 0.74 | 2.95 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | conditional | 1.56 | 0.90 | 0.82 | 2.58 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | conditional | 1.05 | 0.90 | 0.13 | 2.22 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | conditional | 2.08 | 0.90 | 1.61 | 2.60 | 1.00 | Virus_name |
| sigma | fixed | sigma | 1.82 | 0.90 | 1.76 | 1.88 | 1.00 |
Base model with no hierarchical effects:
if (!file.exists("../fit_models/indiv_disease_base.rds")) {
disease_m0 <- brm(Disease_YN ~ Host_order,
data=dat_dis, family=bernoulli(link = "logit"),
prior=c(prior(normal(0, 1), class = Intercept),
prior(normal(0, 1.5), class = b)),
iter=4000, thin=1,
control=list(adapt_delta=0.8, max_treedepth=10), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(disease_m0, "../fit_models/indiv_disease_base.rds")
} else { disease_m0 <- readRDS("../fit_models/indiv_disease_base.rds")}
model <- disease_m0
# summary(model, prob=0.9)
# save memory
# rm(disease_m0)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Mean | CI | CI_low | CI_high | Rhat |
|---|---|---|---|---|---|
| b_Intercept | -0.57 | 0.90 | -0.71 | -0.44 | 1.00 |
| b_Host_orderChiroptera | 0.40 | 0.90 | 0.24 | 0.56 | 1.00 |
p1 <- cond_eff_plot(model) + ylab("Conditional log-odds of disease") +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
Sample size: 2444
With no hierarchical effects and only the bat / rodent binary predictor, the model suggests that bats are more likely than rodents to suffer overt disease.
However, once include hierarchical effects for host and virus species, this effect goes away:
if (!file.exists("../fit_models/indiv_disease_base_phylo.rds")) {
disease_m1 <- brm(Disease_YN ~ Host_order +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV),
data=dat_dis, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.99, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(disease_m1, "../fit_models/indiv_disease_base_phylo.rds")
} else { disease_m1 <- readRDS("../fit_models/indiv_disease_base_phylo.rds")}
model <- disease_m1
# summary(model, prob=0.9)
# save memory
# rm(disease_m1)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.47 | 0.90 | -1.24 | 2.18 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.69 | 0.90 | -2.37 | 1.04 | 1.00 | |
| sd_Host_name__Intercept | random | 3.41 | 0.90 | 2.69 | 4.18 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 1.15 | 0.90 | 0.11 | 2.50 | 1.00 | Host_Upham |
| sd_Virus_ICTV__Intercept | random | 1.95 | 0.90 | 0.76 | 3.12 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 2.58 | 0.90 | 1.80 | 3.41 | 1.00 | Virus_name |
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) +
ylab("Conditional log-odds of disease") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_disease_base_phylo_comboplot.pdf", width=8, height=3)
Sample size: 2444
Full model adjusts for evolutionary predictors, dose, and inoculation route
if (!file.exists("../fit_models/indiv_disease_full.rds")) {
disease_full <- brm(Disease_YN ~ Host_order +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_dis, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.95, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(disease_full, "../fit_models/indiv_disease_full.rds")
} else { disease_full <- readRDS("../fit_models/indiv_disease_full.rds")}
model <- disease_full
# summary(model, prob=0.9)
# save memory
# rm(disease_full)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.18 | 0.90 | -1.64 | 2.01 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.54 | 0.90 | -2.20 | 1.14 | 1.00 | |
| b_ses_pd | fixed | -1.46 | 0.90 | -2.97 | 0.04 | 1.00 | |
| b_EvoIso | fixed | 0.70 | 0.90 | -0.55 | 1.99 | 1.00 | |
| b_Dose_mass | fixed | 0.71 | 0.90 | 0.26 | 1.18 | 1.00 | |
| sd_Dose_unit__Intercept | random | 2.30 | 0.90 | 1.48 | 3.24 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 2.77 | 0.90 | 2.10 | 3.47 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 1.01 | 0.90 | 0.08 | 2.36 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.62 | 0.90 | 0.98 | 2.45 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 3.17 | 0.90 | 2.06 | 4.21 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.19 | 0.90 | 0.14 | 2.39 | 1.00 | Virus_name |
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,15)) +
ylab("Conditional log-odds of disease") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_disease_full_comboplot.pdf", width=8, height=3)
Sample size: 1500
dat_mort <- dat
dat_mort <- dat_mort[dat_mort$Susceptible_YN==1,]
dat_mort <- dat_mort[!is.na(dat_mort$Susceptible_YN),]
dat_mort <- dat_mort[!is.na(dat_mort$Mortality_YN),]
dat_mort$ses_pd <- scale_half(dat_mort$ses_pd)
dat_mort$EvoIso <- scale_half(dat_mort$EvoIso)
dat_mort <- dat_mort %>% group_by(Dose_unit) %>% mutate(Dose_mass = scale_half(log(Dose_mass)))
dat_mort$Host_name <- dat_mort$Host_Upham
dat_mort$Virus_name <- dat_mort$Virus_ICTV
if (!file.exists("../fit_models/indiv_mort_full.rds")) {
mort_full <- brm(Mortality_YN ~ Host_order +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_mort, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1,
cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.95, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(mort_full, "../fit_models/indiv_mort_full.rds")
} else { mort_full <- readRDS("../fit_models/indiv_mort_full.rds")}
model <- mort_full
# summary(model, prob=0.9)
# save memory
# rm(mort_full)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | -1.59 | 0.90 | -3.22 | 0.10 | 1.00 | |
| b_Host_orderChiroptera | fixed | 1.31 | 0.90 | -0.40 | 3.01 | 1.00 | |
| b_ses_pd | fixed | -2.25 | 0.90 | -3.32 | -1.18 | 1.00 | |
| b_EvoIso | fixed | 0.43 | 0.90 | -0.79 | 1.74 | 1.00 | |
| b_Dose_mass | fixed | 5.75e-03 | 0.90 | -0.52 | 0.54 | 1.00 | |
| sd_Dose_unit__Intercept | random | 1.75 | 0.90 | 1.05 | 2.65 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 2.28 | 0.90 | 1.66 | 2.92 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.88 | 0.90 | 0.08 | 2.03 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.92 | 0.90 | 1.15 | 2.87 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 0.59 | 0.90 | 0.04 | 1.50 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.59 | 0.90 | 0.05 | 1.42 | 1.00 | Virus_name |
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-18,18)) +
ylab("Conditional log-odds of mortality") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar2 <- intervals_plot(model)
plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_mortality_full_comboplot.pdf", width=8, height=3)
Sample size: 927
Fitting an ordinal regression (a.k.a. cumulative link) model for disease severity
dat_severity <- dat
dat_severity <- dat_severity[dat_severity$Susceptible_YN==1,]
dat_severity <- dat_severity[!is.na(dat_severity$Susceptible_YN),]
dat_severity <- dat_severity[!is.na(dat_severity$severity_rank),]
dat_severity$ses_pd <- scale_half(dat_severity$ses_pd)
dat_severity$EvoIso <- scale_half(dat_severity$EvoIso)
dat_severity <- dat_severity %>% group_by(Dose_unit) %>% mutate(Dose_mass = scale_half(log(Dose_mass)))
dat_severity$Host_name <- dat_severity$Host_Upham
dat_severity$Virus_name <- dat_severity$Virus_ICTV
if (!file.exists("../fit_models/indiv_severity_full.rds")) {
severity_full <- brm(severity_rank ~ Host_order +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_severity, family=cumulative(probit),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.98, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(severity_full, "../fit_models/indiv_severity_full.rds")
} else { severity_full <- readRDS("../fit_models/indiv_severity_full.rds")}
model <- severity_full
# summary(model, prob=0.9)
# save memory
# rm(severity_full)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept[1] | fixed | -0.52 | 0.90 | -1.64 | 0.64 | 1.00 | |
| b_Intercept[2] | fixed | 0.02 | 0.90 | -1.10 | 1.19 | 1.00 | |
| b_Intercept[3] | fixed | 0.32 | 0.90 | -0.80 | 1.48 | 1.00 | |
| b_Intercept[4] | fixed | 0.75 | 0.90 | -0.37 | 1.92 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.06 | 0.90 | -1.32 | 1.30 | 1.00 | |
| b_ses_pd | fixed | -1.16 | 0.90 | -2.05 | -0.27 | 1.00 | |
| b_EvoIso | fixed | 0.22 | 0.90 | -0.43 | 0.91 | 1.00 | |
| b_Dose_mass | fixed | 0.60 | 0.90 | 0.40 | 0.80 | 1.00 | |
| sd_Dose_unit__Intercept | random | 1.72 | 0.90 | 1.04 | 2.61 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 1.48 | 0.90 | 0.99 | 1.95 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.90 | 0.90 | 0.08 | 1.97 | 1.01 | Host_Upham |
| sd_Route_type__Intercept | random | 1.00 | 0.90 | 0.51 | 1.75 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 2.01 | 0.90 | 1.39 | 2.67 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.40 | 0.90 | 0.03 | 0.98 | 1.00 | Virus_name |
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) +
ylab("Conditional log-odds of severity") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar3 <- intervals_plot(model)
plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_severity_full_comboplot.pdf", width=8, height=3)
Sample size: 1500
## Joint plot
global_ylim <- scale_y_continuous(lim=c(-15,15))
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))
joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)
ggsave("../plots_tables/individual_level_conditionals_intervals_full.pdf", joint_plot, width=10.5, height=5.5)
if (!file.exists("../fit_models/indiv_disease_resMatch.rds")) {
disease_resMatch <- brm(Disease_YN ~ Host_order + Reservoir_match +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_dis, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.95, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(disease_resMatch, "../fit_models/indiv_disease_resMatch.rds")
} else { disease_resMatch <- readRDS("../fit_models/indiv_disease_resMatch.rds")}
model <- disease_resMatch
# summary(model, prob=0.9)
# save memory
# rm(disease_resMatch)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | -0.50 | 0.90 | -2.68 | 1.67 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.54 | 0.90 | -2.19 | 1.13 | 1.00 | |
| b_Reservoir_match | fixed | 1.09 | 0.90 | -0.35 | 2.57 | 1.00 | |
| b_ses_pd | fixed | -1.26 | 0.90 | -2.85 | 0.28 | 1.00 | |
| b_EvoIso | fixed | 0.58 | 0.90 | -0.77 | 1.87 | 1.00 | |
| b_Dose_mass | fixed | 0.67 | 0.90 | 0.21 | 1.12 | 1.00 | |
| sd_Dose_unit__Intercept | random | 2.18 | 0.90 | 1.38 | 3.11 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 2.60 | 0.90 | 1.93 | 3.33 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.90 | 0.90 | 0.07 | 2.19 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.67 | 0.90 | 1.01 | 2.50 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 2.59 | 0.90 | 0.87 | 3.97 | 1.01 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.49 | 0.90 | 0.20 | 2.76 | 1.01 | Virus_name |
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,15)) +
ylab("Conditional log-odds of disease") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_disease_resMatch_comboplot.pdf", width=8, height=3)
Sample size: 1391
if (!file.exists("../fit_models/indiv_mortality_resMatch.rds")) {
mortality_resMatch <- brm(Mortality_YN ~ Host_order + Reservoir_match +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_mort, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.95, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(mortality_resMatch, "../fit_models/indiv_mortality_resMatch.rds")
} else { mortality_resMatch <- readRDS("../fit_models/indiv_mortality_resMatch.rds")}
model <- mortality_resMatch
# summary(model, prob=0.9)
# save memory
# rm(mortality_resMatch)
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | -1.62 | 0.90 | -3.45 | 0.20 | 1.00 | |
| b_Host_orderChiroptera | fixed | 1.03 | 0.90 | -0.59 | 2.64 | 1.00 | |
| b_Reservoir_match | fixed | 0.29 | 0.90 | -1.06 | 1.69 | 1.00 | |
| b_ses_pd | fixed | -1.97 | 0.90 | -3.15 | -0.71 | 1.00 | |
| b_EvoIso | fixed | 0.51 | 0.90 | -0.80 | 1.90 | 1.00 | |
| b_Dose_mass | fixed | 0.02 | 0.90 | -0.52 | 0.57 | 1.00 | |
| sd_Dose_unit__Intercept | random | 1.70 | 0.90 | 1.00 | 2.56 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 2.19 | 0.90 | 1.60 | 2.82 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.75 | 0.90 | 0.06 | 1.81 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.97 | 0.90 | 1.18 | 2.95 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 0.66 | 0.90 | 0.05 | 1.66 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.74 | 0.90 | 0.07 | 1.69 | 1.00 | Virus_name |
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Mortality_YN, pp[1:500, ])
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,15)) +
ylab("Conditional log-odds of mortality") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar2 <- intervals_plot(model)
plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_mortality_resMatch_comboplot.pdf", width=8, height=3)
Sample size: 851
if (!file.exists("../fit_models/indiv_severity_resMatch.rds")) {
severity_resMatch <- brm(severity_rank ~ Host_order + Reservoir_match +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_severity, family=cumulative(probit),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.98, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(severity_resMatch, "../fit_models/indiv_severity_resMatch.rds")
} else { severity_resMatch <- readRDS("../fit_models/indiv_severity_resMatch.rds")}
model <- severity_resMatch
# summary(model, prob=0.9)
# save memory
# rm(severity_resMatch)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$severity_rank, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept[1] | fixed | 0.01 | 0.90 | -1.37 | 1.40 | 1.00 | |
| b_Intercept[2] | fixed | 0.50 | 0.90 | -0.89 | 1.89 | 1.00 | |
| b_Intercept[3] | fixed | 0.80 | 0.90 | -0.59 | 2.20 | 1.00 | |
| b_Intercept[4] | fixed | 1.24 | 0.90 | -0.14 | 2.63 | 1.00 | |
| b_Host_orderChiroptera | fixed | -8.57e-03 | 0.90 | -1.22 | 1.27 | 1.00 | |
| b_Reservoir_match | fixed | 0.72 | 0.90 | -0.19 | 1.66 | 1.00 | |
| b_ses_pd | fixed | -1.02 | 0.90 | -2.08 | -0.04 | 1.00 | |
| b_EvoIso | fixed | 0.16 | 0.90 | -0.64 | 0.98 | 1.00 | |
| b_Dose_mass | fixed | 0.59 | 0.90 | 0.39 | 0.79 | 1.00 | |
| sd_Dose_unit__Intercept | random | 1.69 | 0.90 | 1.01 | 2.56 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 1.43 | 0.90 | 0.98 | 1.88 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.72 | 0.90 | 0.05 | 1.81 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.00 | 0.90 | 0.51 | 1.76 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 1.78 | 0.90 | 0.68 | 2.71 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.74 | 0.90 | 0.09 | 1.54 | 1.00 | Virus_name |
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,15)) +
ylab("Conditional log-odds of severity") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar3 <- intervals_plot(model)
plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_severity_resMatch_comboplot.pdf", width=8, height=3)
Sample size: 1391
## Joint plot
global_ylim <- scale_y_continuous(lim=c(-18,18))
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))
joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)
ggsave("../plots_tables/individual_level_conditionals_intervals_resMatch.pdf", joint_plot, width=10.5, height=5.5)
lyssaviruses <- viruses$Virus_ICTV[viruses$genus=="Lyssavirus"] %>% unique()
dat_dis_nolys <- dat_dis[!dat_dis$Virus_ICTV%in%lyssaviruses,]
if (!file.exists("../fit_models/indiv_disease_nolys.rds")) {
disease_nolys <- brm(Disease_YN ~ Host_order +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_dis_nolys, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.9, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(disease_nolys, "../fit_models/indiv_disease_nolys.rds")
} else { disease_nolys <- readRDS("../fit_models/indiv_disease_nolys.rds")}
model <- disease_nolys
# summary(model, prob=0.9)
# save memory
# rm(disease_nolys)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | 0.23 | 0.90 | -1.61 | 2.06 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.62 | 0.90 | -2.34 | 1.20 | 1.00 | |
| b_ses_pd | fixed | -1.15 | 0.90 | -3.16 | 0.82 | 1.00 | |
| b_EvoIso | fixed | 0.77 | 0.90 | -0.73 | 2.33 | 1.00 | |
| b_Dose_mass | fixed | 1.05 | 0.90 | 0.37 | 1.76 | 1.00 | |
| sd_Dose_unit__Intercept | random | 1.46 | 0.90 | 0.31 | 2.65 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 3.26 | 0.90 | 2.53 | 4.05 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 1.04 | 0.90 | 0.08 | 2.49 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.24 | 0.90 | 0.50 | 2.18 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 1.77 | 0.90 | 0.21 | 3.43 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 2.87 | 0.90 | 1.69 | 3.97 | 1.00 | Virus_name |
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-18,18)) +
ylab("Conditional log-odds of disease") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_disease_nolys_comboplot.pdf", width=8, height=3)
Sample size: 1425
lyssaviruses <- viruses$Virus_ICTV[viruses$genus=="Lyssavirus"] %>% unique()
dat_mort_nolys <- dat_mort[!dat_mort$Virus_ICTV%in%lyssaviruses,]
if (!file.exists("../fit_models/indiv_mort_nolys.rds")) {
mort_nolys <- brm(Mortality_YN ~ Host_order +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_mort_nolys, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(mort_nolys, "../fit_models/indiv_mort_nolys.rds")
} else { mort_nolys <- readRDS("../fit_models/indiv_mort_nolys.rds")}
model <- mort_nolys
# summary(model, prob=0.9)
# save memory
# rm(mort_nolys)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | -0.91 | 0.90 | -2.56 | 0.78 | 1.00 | |
| b_Host_orderChiroptera | fixed | 0.19 | 0.90 | -1.67 | 2.09 | 1.00 | |
| b_ses_pd | fixed | -1.55 | 0.90 | -3.16 | 0.13 | 1.00 | |
| b_EvoIso | fixed | 1.40 | 0.90 | -0.25 | 3.11 | 1.00 | |
| b_Dose_mass | fixed | 0.28 | 0.90 | -0.60 | 1.18 | 1.00 | |
| sd_Dose_unit__Intercept | random | 1.32 | 0.90 | 0.18 | 2.55 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 2.18 | 0.90 | 1.45 | 2.95 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 1.05 | 0.90 | 0.09 | 2.37 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.56 | 0.90 | 0.65 | 2.67 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 0.84 | 0.90 | 0.07 | 2.00 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 0.96 | 0.90 | 0.10 | 2.02 | 1.00 | Virus_name |
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-18,10)) +
ylab("Conditional log-odds of mortality") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar2 <- intervals_plot(model)
plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_mortality_nolys_comboplot.pdf", width=8, height=3)
Sample size: 742
lyssaviruses <- viruses$Virus_ICTV[viruses$genus=="Lyssavirus"] %>% unique()
dat_severity_nolys <- dat_severity[!dat_severity$Virus_ICTV%in%lyssaviruses,]
if (!file.exists("../fit_models/indiv_severity_nolys.rds")) {
severity_nolys <- brm(severity_rank ~ Host_order +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_severity_nolys, family=cumulative(probit),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(severity_nolys, "../fit_models/indiv_severity_nolys.rds")
} else { severity_nolys <- readRDS("../fit_models/indiv_severity_nolys.rds")}
model <- severity_nolys
# summary(model, prob=0.9)
# save memory
# rm(severity_nolys)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept[1] | fixed | -1.17 | 0.90 | -2.45 | 0.11 | 1.00 | |
| b_Intercept[2] | fixed | -0.34 | 0.90 | -1.62 | 0.94 | 1.00 | |
| b_Intercept[3] | fixed | 0.08 | 0.90 | -1.20 | 1.37 | 1.00 | |
| b_Intercept[4] | fixed | 0.65 | 0.90 | -0.63 | 1.93 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.28 | 0.90 | -1.75 | 1.32 | 1.00 | |
| b_ses_pd | fixed | -1.23 | 0.90 | -2.69 | 0.25 | 1.00 | |
| b_EvoIso | fixed | 0.50 | 0.90 | -0.42 | 1.41 | 1.00 | |
| b_Dose_mass | fixed | 0.42 | 0.90 | 0.16 | 0.69 | 1.00 | |
| sd_Dose_unit__Intercept | random | 1.28 | 0.90 | 0.38 | 2.24 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 1.98 | 0.90 | 1.38 | 2.58 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 1.08 | 0.90 | 0.09 | 2.42 | 1.01 | Host_Upham |
| sd_Route_type__Intercept | random | 1.35 | 0.90 | 0.65 | 2.23 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 1.74 | 0.90 | 0.43 | 2.82 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.11 | 0.90 | 0.14 | 2.10 | 1.00 | Virus_name |
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,10)) +
ylab("Conditional log-odds of severity") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar3 <- intervals_plot(model)
plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_severity_nolys_comboplot.pdf", width=8, height=3)
Sample size: 1425
## Joint plot
global_ylim <- scale_y_continuous(lim=c(-18,10))
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))
joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)
ggsave("../plots_tables/individual_level_conditionals_intervals_nolys.pdf", joint_plot, width=10.5, height=5.5)
dat_dis_heteros <- dat_dis[dat_dis$Reservoir_match%in%0,]
if (!file.exists("../fit_models/indiv_disease_heteroHosts.rds")) {
disease_heteros <- brm(Disease_YN ~ Host_order +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_dis_heteros, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.95, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(disease_heteros, "../fit_models/indiv_disease_heteroHosts.rds")
} else { disease_heteros <- readRDS("../fit_models/indiv_disease_heteroHosts.rds")}
model <- disease_heteros
# summary(model, prob=0.9)
# save memory
# rm(disease_heteros)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | -0.17 | 0.90 | -2.12 | 1.78 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.67 | 0.90 | -2.42 | 1.08 | 1.00 | |
| b_ses_pd | fixed | -0.62 | 0.90 | -2.78 | 1.51 | 1.00 | |
| b_EvoIso | fixed | 1.16 | 0.90 | -0.39 | 2.89 | 1.00 | |
| b_Dose_mass | fixed | 1.02 | 0.90 | 0.30 | 1.78 | 1.00 | |
| sd_Dose_unit__Intercept | random | 1.17 | 0.90 | 0.17 | 2.31 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 2.66 | 0.90 | 1.94 | 3.45 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.78 | 0.90 | 0.06 | 1.93 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.60 | 0.90 | 0.84 | 2.55 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 1.05 | 0.90 | 0.11 | 2.26 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.25 | 0.90 | 0.16 | 2.44 | 1.00 | Virus_name |
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,12)) +
ylab("Conditional log-odds of disease") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) +
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_disease_heteroHosts_comboplot.pdf", width=8, height=3)
Sample size: 872
dat_mort_heteros <- dat_mort[dat_mort$Reservoir_match%in%0,]
if (!file.exists("../fit_models/indiv_mort_heteroHosts.rds")) {
mort_heteros <- brm(Mortality_YN ~ Host_order +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_mort_heteros, family=bernoulli(link = "logit"),
prior=custom_prior,
iter=5000, thin=1,
cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.95, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(mort_heteros, "../fit_models/indiv_mort_heteroHosts.rds")
} else { mort_heteros <- readRDS("../fit_models/indiv_mort_heteroHosts.rds")}
model <- mort_heteros
# summary(model, prob=0.9)
# save memory
# rm(mort_heteros)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept | fixed | -1.06 | 0.90 | -2.83 | 0.72 | 1.00 | |
| b_Host_orderChiroptera | fixed | 0.11 | 0.90 | -1.71 | 1.91 | 1.00 | |
| b_ses_pd | fixed | 0.18 | 0.90 | -1.98 | 2.30 | 1.00 | |
| b_EvoIso | fixed | 1.06 | 0.90 | -0.79 | 3.07 | 1.00 | |
| b_Dose_mass | fixed | 0.37 | 0.90 | -0.61 | 1.36 | 1.00 | |
| sd_Dose_unit__Intercept | random | 1.05 | 0.90 | 0.10 | 2.34 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 1.45 | 0.90 | 0.51 | 2.40 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.92 | 0.90 | 0.09 | 2.10 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.78 | 0.90 | 0.93 | 2.80 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 0.85 | 0.90 | 0.07 | 2.05 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.45 | 0.90 | 0.31 | 2.57 | 1.00 | Virus_name |
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-18,12)) +
ylab("Conditional log-odds of mortality") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar2 <- intervals_plot(model)
plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_mortality_heteroHosts_comboplot.pdf", width=8, height=3)
Sample size: 358
dat_severity_heteros <- dat_severity[dat_severity$Reservoir_match%in%0,]
if (!file.exists("../fit_models/indiv_severity_heteroHosts.rds")) {
severity_heteros <- brm(severity_rank ~ Host_order +
ses_pd + EvoIso +
Dose_mass + (1|Dose_unit) +
(1|Host_name) + (1|Virus_name) +
(1|Host_Upham) + (1|Virus_ICTV) +
(1|Route_type),
data=dat_severity_heteros, family=cumulative(probit),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
control=list(adapt_delta=0.98, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(severity_heteros, "../fit_models/indiv_severity_heteroHosts.rds")
} else { severity_heteros <- readRDS("../fit_models/indiv_severity_heteroHosts.rds")}
model <- severity_heteros
# summary(model, prob=0.9)
# save memory
# rm(severity_heteros)
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])
param_tab(model)
| Parameter | Effects | Mean | CI | CI_low | CI_high | Rhat | Group |
|---|---|---|---|---|---|---|---|
| b_Intercept[1] | fixed | -0.47 | 0.90 | -1.80 | 0.87 | 1.00 | |
| b_Intercept[2] | fixed | -0.27 | 0.90 | -1.61 | 1.06 | 1.00 | |
| b_Intercept[3] | fixed | 0.08 | 0.90 | -1.26 | 1.41 | 1.00 | |
| b_Intercept[4] | fixed | 0.65 | 0.90 | -0.68 | 1.98 | 1.00 | |
| b_Host_orderChiroptera | fixed | -0.39 | 0.90 | -1.73 | 0.98 | 1.00 | |
| b_ses_pd | fixed | -0.60 | 0.90 | -2.45 | 1.21 | 1.00 | |
| b_EvoIso | fixed | 0.61 | 0.90 | -0.64 | 2.02 | 1.00 | |
| b_Dose_mass | fixed | 0.47 | 0.90 | 0.16 | 0.78 | 1.00 | |
| sd_Dose_unit__Intercept | random | 0.94 | 0.90 | 0.13 | 1.89 | 1.00 | Dose_unit |
| sd_Host_name__Intercept | random | 1.30 | 0.90 | 0.75 | 1.85 | 1.00 | Host_name |
| sd_Host_Upham__Intercept | random | 0.75 | 0.90 | 0.07 | 1.74 | 1.00 | Host_Upham |
| sd_Route_type__Intercept | random | 1.57 | 0.90 | 0.88 | 2.46 | 1.00 | Route_type |
| sd_Virus_ICTV__Intercept | random | 0.93 | 0.90 | 0.10 | 2.03 | 1.00 | Virus_ICTV |
| sd_Virus_name__Intercept | random | 1.17 | 0.90 | 0.31 | 2.04 | 1.00 | Virus_name |
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-10,5)) +
ylab("Conditional log-odds of severity") + xlab("") +
theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar3 <- intervals_plot(model)
plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))
# ggsave("../plots_tables/individual_level_severity_heteroHosts_comboplot.pdf", width=8, height=3)
Sample size: 872
## Joint plot
global_ylim <- scale_y_continuous(lim=c(-15,10))
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))
joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)
ggsave("../plots_tables/individual_level_conditionals_intervals_heteroHosts.pdf", joint_plot, width=10.5, height=5.5)
s_dis <- readRDS("../fit_models/spec_disease_full.rds")
s_mort <- readRDS("../fit_models/spec_mort_full.rds")
s_sev <- readRDS("../fit_models/spec_severity_full.rds")
i_dis <- readRDS("../fit_models/indiv_disease_full.rds")
i_mort <- readRDS("../fit_models/indiv_mort_full.rds")
i_sev <- readRDS("../fit_models/indiv_severity_full.rds")
padding <- theme(plot.margin = unit(c(t=0.1, r=0.1, b=0, l=0), "cm"))
cp <- egg::ggarrange(
(compare_posteriors(s_dis, i_dis) + padding + theme(axis.text.y=element_text(size=12)) + ggtitle("A) Disease presence")),
(compare_posteriors(s_mort, i_mort) + padding + theme(axis.text.y=element_blank()) + ggtitle("B) Mortality")),
(compare_posteriors(s_sev, i_sev) + padding + theme(axis.text.y=element_blank(), legend.position = "right") + ggtitle("C) Severity")),
ncol=3, padding=0.5)
ggsave("../plots_tables/combined_species_individual_intervals.pdf", cp, width=11.5, height=3.5)
| Host_order | mean_cfr | median_cfr | min_cfr | max_cfr |
|---|---|---|---|---|
| Rodentia | 7.57 | 1.17 | 0 | 100 |
| Chiroptera | 62.90 | 80.92 | 0 | 100 |
Comparing CFR for all viruses, regardless of host origin, and excluding bat viruses.
dat_cfr$Virus_name <- dat_cfr$Virus_ICTV
dat_cfr$Host_name <- dat_cfr$Host_Upham
# do bats and rodents differ in mean human CFR of inoculated viruses?
# model with no taxonomic / phylogenetic effects - rodents infected with lower human cfr viruses
if (!file.exists("../fit_models/cfr_base1.rds")) {
cfr_base1 <- brm(CFR_avg/100 ~ Host_order+
(1|PaperID),
data=dat_cfr, family=zero_one_inflated_beta(),
prior=custom_prior,
iter=4000, thin=1,
control=list(adapt_delta=0.85, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(cfr_base1, "../fit_models/cfr_base1.rds")
} else { cfr_base1 <- readRDS("../fit_models/cfr_base1.rds")}
# cfr_base1
# pp_check(cfr_base1)
# cond_eff_plot(cfr_base1)
# param_tab(cfr_base1)
posterior <- as.matrix(cfr_base1)
color_scheme_set("teal")
ar1 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean",
point_size = 3.5, outer_size=1, inner_size=3)
# ar1
# model with only viral taxonomic hierarcical effects
if (!file.exists("../fit_models/cfr_base2.rds")) {
cfr_base2 <- brm(CFR_avg/100 ~ Host_order +
(1|Virus_name) +
(1|Virus_ICTV) +
(1|PaperID),
data=dat_cfr, family=zero_one_inflated_beta(),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(cfr_base2, "../fit_models/cfr_base2.rds")
} else { cfr_base2 <- readRDS("../fit_models/cfr_base2.rds")}
# cfr_base2
# pp_check(cfr_base2, ndraws = 200)
# cond_eff_plot(cfr_base2)
# param_tab(cfr_base2)
color_scheme_set("teal")
posterior <- as.matrix(cfr_base2)
ar2 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean",
point_size = 3.5, outer_size=1, inner_size=3)
# ar2
# model with no viral effects, but removing bat-reservoired viruses
bat_res_viruses <- vtraits$Virus_ICTV[vtraits$Reservoir%in%"Chiroptera"] %>% unique()
dat_cfr_nobatv <- dat_cfr[!dat_cfr$Virus_ICTV%in%bat_res_viruses,]
if (!file.exists("../fit_models/cfr_base3.rds")) {
cfr_nobatv <- brm(CFR_avg/100 ~ Host_order +
(1|PaperID),
data=dat_cfr_nobatv, family=zero_one_inflated_beta(),
prior=custom_prior,
iter=4000, thin=1,
control=list(adapt_delta=0.85, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(cfr_nobatv, "../fit_models/cfr_base3.rds")
} else { cfr_nobatv <- readRDS("../fit_models/cfr_base3.rds")}
# cfr_nobatv
# pp_check(cfr_nobatv, ndraws=500)
# cond_eff_plot(cfr_nobatv)
# param_tab(cfr_nobatv)
posterior <- as.matrix(cfr_nobatv)
ar3 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean",
point_size = 3.5, outer_size=1, inner_size=3)
# ar3
# # combo plot
# p1 <- plot_grid( ar1 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nAll viruses"),
# ar3 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nExcluding bat-reservoired viruses"),
# ar2 + xlim(-3, 2) + scale_y_discrete(labels = "With virus effects\nAll viruses"),
# ncol=1, align = "v")
# p2 <- add_sub(p1, "Host order (Chiroptera) effect", x = 0, hjust = -3.5, size=10)
# ggdraw(p2)
# ggsave("../plots_tables/CFR_models_intervals.pdf", width=6, height=4)
### Models with only heterologous hosts
dat_cfr_het <- left_join(dat_cfr, dat %>% select(Host_Upham, Virus_ICTV, Reservoir_match, PaperID) %>% unique())
dat_cfr_het <- dat_cfr_het[dat_cfr_het$Reservoir_match==0,]
# model with no virus hierarchical effects
if (!file.exists("../fit_models/cfr_het1.rds")) {
cfr_het1 <- brm(CFR_avg/100 ~ Host_order + (1|PaperID),
data=dat_cfr_het, family=zero_one_inflated_beta(),
prior=custom_prior,
iter=4000, thin=1,
control=list(adapt_delta=0.85, max_treedepth=12), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(cfr_het1, "../fit_models/cfr_het1.rds")
} else { cfr_het1 <- readRDS("../fit_models/cfr_het1.rds")}
# cfr_het1
# pp_check(cfr_het1, ndraws=300)
# cond_eff_plot(cfr_het1)
# param_tab(cfr_het1)
posterior <- as.matrix(cfr_het1)
ar4 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean",
point_size = 3.5, outer_size=1, inner_size=3)
# ar4
# model with only viral taxonomic hierarcical effects
if (!file.exists("../fit_models/cfr_het2.rds")) {
cfr_het2 <- brm(CFR_avg/100 ~ Host_order +
(1|Virus_name) +
(1|Virus_ICTV) +
(1|PaperID),
data=dat_cfr_het, family=zero_one_inflated_beta(),
prior=custom_prior,
iter=5000, thin=1, cov_ranef = list(Virus_ICTV=virus_cov),
control=list(adapt_delta=0.90, max_treedepth=16), cores=4,
save_pars=save_pars(all=TRUE))
saveRDS(cfr_het2, "../fit_models/cfr_het2.rds")
# 0 div trans
} else { cfr_het2 <- readRDS("../fit_models/cfr_het2.rds")}
# pairs(cfr_het2)
# cond_eff_plot(cfr_het2)
# param_tab(cfr_het2)
color_scheme_set("teal")
posterior <- as.matrix(cfr_het2)
ar5 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean",
point_size = 3.5, outer_size=1, inner_size=3)
# ar5
# combo plot
p1 <- plot_grid( ar1 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nAll viruses"),
ar3 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nExcluding bat-reservoired viruses"),
ar2 + xlim(-3, 2) + scale_y_discrete(labels = "With virus effects\nAll viruses"),
ar4 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nHeterologous hosts"),
ar5 + xlim(-3, 2) + scale_y_discrete(labels = "With virus effects\nHeterologous hosts"),
ncol=1, align = "v")
p2 <- add_sub(p1, "Host order (Chiroptera) effect", x = 0, hjust = -3.5, size=10)
ggdraw(p2)
ggsave("../plots_tables/CFR_models_intervals.pdf", width=6, height=6)
# dat_cfr %>% ungroup() %>% select(Virus_ICTV, CFR_avg) %>% unique() %>% arrange(CFR_avg) %>% View()